The alternating-direction solver algorithm was first proposed by Jim Douglas in 1956 in the context of non-stationary simulations using the finite difference method [1].
The stationary version for the isogeometric finite element method was proposed by Longfei Gao in his doctoral dissertation under the supervision of prof. Victor Manuel Calo [2].
Consider the system of equations we got in the first chapter (related to the bitmap projection problem). As a result of the generation of the system of equations, we got a matrix written in the form of a Kronecker product of two five-diagonal matrices, and the right-hand side counted for individual testing B-splines
\( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{20} & \frac{13}{120} & \frac{1}{120} \\ \end{bmatrix} \otimes \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \begin{bmatrix} u_{1,1} \\ u_{2,1} \\ u_{3,1} \\ \vdots \\ u_{k,l} \\ \vdots \\ u_{N_{x-2},N_y} \\ u_{N_{x-1},N_y} \\ u_{N_x,N_y} \\ \end{bmatrix} \\ = \begin{bmatrix} \int BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_1(y)} dxdy \\ \int BITMAP(x,y) {\color{blue}B^x_2(x)*B^y_1(y)} dxdy \\ \int BITMAP(x,y) {\color{blue}B^x_3(x)*B^y_1(y)} dxdy \\ \vdots \\ \int BITMAP(x,y) {\color{blue}B^x_k(x)*B^y_l(y)} dxdy \\ \vdots \\ \int BITMAP(x,y) {\color{blue}B^x_{N_x-2}(x)*B^y_{N_y}(y)} dxdy \\ \int BITMAP(x,y) {\color{blue}B^x_{N_x-1}(x)*B^y_{N_y}(y)} dxdy \\ \int BITMAP(x,y) {\color{blue}B^x_{N_x}(x)*B^y_{N_y}(y)} dxdy \\ \end{bmatrix} \)
The solution of a system of equations in which the matrix has the structure of a Kronecker product is possible in a very short time. What does this mean in a very short time? A computational cost is expressed in the number of operations, such as the multiplication or addition of numbers necessary to solve a system of equations. In the case of a system of equations in which the matrix has the structure of a Kronecker product, it is possible to solve the system of equations using an algorithm in which the number of operations equals \( const*N \) where N is the number of unknowns (the number of bitmap approximation coefficients on the mesh, expressed exactly by \( N=N_x*N_y \) where \( N_x,N_y \) is the mesh dimension, and \( const \) stands for some fixed number.
The algorithm used in this case is called a alternating directions solver algorithm. We consider two steps in the solution process. The first step is to take the first of the sub-matrix and arrange the vector of unknowns and the vector of right-hand sides into many sub-vectors, one vector for each column of computational mesh elements. It is illustrated in Fig. 1, and in the formula below.
\( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \begin{bmatrix} w_{1,1} & w_{1,2} & \cdots & w_{1,N_{y-1}} & w_{1,N_y} \\ w_{2,1} & w_{2,2} & \cdots & w_{2,N_{y-1}} & w_{2,N_y} \\ w_{3,1} & w_{3,2} & \cdots & w_{3,N_{y-1}} & w_{3,N_y} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ w_{N_{x-2},1} & w_{N_{x-2},2} & \cdots & w_{N_{x-2},N_{y-1}} & w_{N_{x-2},N_y} \\ w_{N_{x-1},1} & w_{N_{x-1},2} & \cdots & w_{N_{x-1},N_{y-1}} & w_{N_{x-1},N_y} \\ w_{N_x,1} & w_{N_x,2} & \cdots & w_{N_x,N_{y-1}} & w_{N_x,N_y} \\ \end{bmatrix} \\ = \begin{bmatrix} \int BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_1(y)} dxdy & \cdots & \int BITMAP(x,y) {\color{blue}B^x_1(x)*B^y_{N_y}(y)}dxdy \\ \int BITMAP(x,y) {\color{blue}B^x_2(x)*B^y_1(y)} dxdy & \cdots & \int BITMAP(x,y) {\color{blue}B^x_2(x)*B^y_{N_y}(y)}dxdy \\ \int BITMAP(x,y) {\color{blue}B^x_3(x)*B^y_1(y)} dxdy & \cdots & \int BITMAP(x,y) {\color{blue}B^x_3(x)*B^y_{N_y}(y)}dxdy \\ \vdots \\ \int BITMAP(x,y) {\color{blue}B^x_{N_x-2}(x)*B^y_1(y)} dxdy & \cdots & \int BITMAP(x,y) {\color{blue}B^x_{N_x-2}(x)*B^y_{N_y}(y)}dxdy \\ \int BITMAP(x,y) {\color{blue}B^x_{N_x-1}(x)*B^y_1(y)} dxdy & \cdots & \int BITMAP(x,y) {\color{blue}B^x_{N_x-1}(x)*B^y_{N_y}(y)}dxdy \\ \int BITMAP(x,y) {\color{blue}B^x_{N_x}(x)*B^y_1(y)} dxdy & \cdots & \int BITMAP(x,y) {\color{blue}B^x_{N_x}(x)*B^y_{N_y}(y)}dxdy \\ \end{bmatrix} \)
We have introduced the auxiliary unknowns here \( w* \) which are used to solve the first system of equations. Original unknowns \( u* \) we will get after solving the second system of equations in which the right side will be auxiliary unknowns \( w* \). So we got a system of equations with a five-diagonal matrix with many right sides. Each right-hand side, corresponds to one column on the element grid, so it has a fixed coordinate \( y \), and the coordinate \( x \) varying from 1 to \( N_x \). The u* unknowns, in which the second indexes change with lines, for example, are similarly ordered \( w_{1,1}, w_{1,2}, ..., w_{1,N_y} \) while the first indexes in the columns change.
When we solve this system of equations, we get solutions \( w* \) and we go to the second step of the variable-direction solver algorithm.
The second step, then, is to take an analogous second Kronecker product matrix, take solutions \( w* \) from the first system of equations, arranging them according to the rows of mesh elements (cf. Fig. 1 ), arranging the unknown unknowns in a similar way \( u* \) solving the obtained system of equations with many right-hand sides.

\( \begin{bmatrix} \frac{1}{20} & \frac{13}{120} & \frac{1}{120} & \cdots \\ \frac{13}{120} & \frac{1}{2} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots \\ \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} & \cdots \\ & \vdots & \vdots & \vdots & \vdots & \vdots & \\ & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{11}{20} & \frac{13}{60} & \frac{1}{120} \\ & & \cdots & \frac{1}{120} & \frac{13}{60} & \frac{1}{2} & \frac{13}{120} \\ & & & \cdots & \frac{1}{120} & \frac{13}{120} & \frac{1}{20} \\ \end{bmatrix} \begin{bmatrix} u_{1,1} & u_{2,1} & \cdots & u_{N_{x-1},1} & u_{N_x,1} \\ u_{1,2} & u_{2,2} & \cdots & u_{N_{x-1},2} & u_{N_x,2} \\ u_{1,3} & u_{2,3} & \cdots & u_{N_{x-1},3} & u_{N_x,3} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ u_{1,N_{y-2}} & u_{2,N_{y-2}} & \cdots & u_{N_{x-1},N_{y-2}} & u_{N_{x},N_{y-2}} \\ u_{1,N_{y-1}} & u_{2,N_{y-1}} & \cdots & u_{N_{x-1},N_{y-1}} & u_{N_{x},N_{y-1}} \\ u_{1,N_y} & u_{2,N_y} & \cdots & u_{N_{x-1},N_{y}} & u_{N_x,N_y} \\ \end{bmatrix} \\ = \begin{bmatrix} w_{1,1} & w_{2,1} & \cdots & w_{N_{x-1},1} & w_{N_x,1} \\ w_{1,2} & w_{2,2} & \cdots & w_{N_{x-1},2} & w_{N_x,2} \\ w_{1,3} & w_{2,3} & \cdots & w_{N_{x-1},3} & w_{N_x,3} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ w_{1,N_{y-2}} & w_{2,N_{y-2}} & \cdots & w_{N_{x-1},N_{y-2}} & w_{N_{x},N_{y-2}} \\ w_{1,N_{y-1}} & w_{2,N_{y-1}} & \cdots & w_{N_{x-1},N_{y-1}} & w_{N_{x},N_{y-1}} \\ w_{1,N_y} & w_{2,N_y} & \cdots & w_{N_{x-1},N_{y}} & w_{N_x,N_y} \\ \end{bmatrix} \)
In this second system of equations, each sub-sector, each right-hand side, corresponds to one row in the element grid, so it has a fixed coordinate x, and a coordinate \( x \), and a coordinate \( y \) varying from 1 to \( N_y \).
The unknowns are similarly ordered \( u* \), in which the first indexes change, for example \( w_{1,1}, w_{2,1}, ..., w_{N_x,1} \) while the second indexes change in the columns.
Each of the above equations can now be solved using the Gaussian elimination algorithm for the band matrix.